//  surl.cc - Solve Unit Resitance Lattice
//
//  Usage:  surl  [-verbose]  [inputfile,default=stdin]
//          surl  -help
//
//  Author: Fred Curtis (fred@f2.org)
//  Written: 2000-01-19 - Lattices only
//  Modified 2000-01-27 - Added hypercubes, simplexes and cross polytopes
//  Modified 2000-01-28 - Added trigrid and hexgrid
//  For some notes on resistance problems see: http://f2.org/maths/resnet/
//
//  Written in C++; compiles with g++ under Solaris and WinNT at least;
//  probably won't compile with Sun's CC since that compiler is broken
//  w.r.t. the scope of variables declare within for loops 
//
//  Note: 2004-11-16 - Compiles with gcc 2.96 under Redhat Linux 7.3 with:
//    gcc -O -lstdc++ -o surl surl.cc
//  
#include <stdio.h>    // for fread()
#include <stdlib.h>   // for atoi()
#include <string.h>
#include <ctype.h>
#include <iostream.h>

static char *progName;
static char *inFile;
static int  lineNo = 1;
static char lineBuf[1024];

#define MAXDIM 30

static void fmtError(char *msg, bool give_help = true) {
  static bool gave_help = false;
  if (msg != 0) {
    cerr << progName << " - Error on " << inFile
         << " line " << lineNo << ": " << msg << "\n"
         << "\t" << lineBuf << "\n";
  }
  if (give_help && !gave_help) {
    gave_help = true;
    cerr <<
"Input lines should be either:\n"
"   Empty lines      - only space characters on the line\n"
"   Comment lines    - first non-blank character is a '#'\n"
"   Data lines       - single line with the following format:\n"
"      SIMPLEX   dim            --- 'dim'-dimensional simplex\n"
"      HYPERCUBE dim  OPP|ADJ   --- 'dim'-dimensional hypercube\n"
"      CROSSPOLY dim  OPP|ADJ   --- 'dim'-dimensional cross polytope\n"
"      LATTICE   min-point  ;  max-point  ;  test-point-1  ;  test-point-2\n"
"                               --- lattice, arbitrary dimensions\n"
"      TRIGRID   min-point  ;  max-point  ;  test-point-1  ;  test-point-2\n"
"                               --- same as 2D lattice plus -ve slope diagonals\n"
"      HEXGRID   min-point  ;  max-point  ;  test-point-1  ;  test-point-2\n"
"                               --- same as 2D lattice less edges to left\n"
"                                   of points where (x+y)=0 mod 2\n"
"Example input lines:\n"
"      ## Resistance between opposite corners of a cube\n"
"      LATTICE  0,0,0 ; 1,1,1 ; 0,0,0 ; 1,1,1\n"
"      HYPERCUBE 3 OPP\n"
"      ## Resistance along an edge of a cube\n"
"      LATTICE  0,0,0 ; 1,1,1 ; 0,0,0 ; 1,0,0\n"
"      HYPERCUBE 3 ADJ\n"
"      ## Resistance between opposite corners of a 4D hyper-cube\n"
"      LATTICE  0,0,0,0 ; 1,1,1,1 ; 0,0,0,0 ; 1,1,1,1\n"
"      HYPERCUBE 4 OPP\n"
"      ## Resistance between the central nodes of a 4x3 lattice\n"
"      LATTICE  -1,-1 ; 2,1 ; 0,0 ; 1,0\n"
"      ## Resistance between the central nodes of a 10x9 lattice\n"
"      LATTICE  -4,-4 ; 5,4 ; 0,0 ; 1,0\n"
    ;
  }
}

static void chkDim(int dim) {
  if (dim > MAXDIM) {
    cerr << progName
	 << ": This feeble program is limited to "
	 << MAXDIM << "-dimensions or less\n";
      fmtError(0, false);
      exit(1);
  }
}

static void usageError() {
  cerr << "Usage: "
       << progName << "  [-verbose]  [inputfile,default=stdin]\n";
  fmtError(0);
  exit(1);
}

static double choose(int n, int r) {
  static double arr[MAXDIM+1][MAXDIM+1];
  if (n < 0 || n > MAXDIM || r < 0 || r > n) {
    cerr << "internal error - choose() called with args n="
         << n << ", r=" << r << "\n";
    exit(1);
  }
  if (arr[n][r] == 0) {
    for (int i = 0; i <= MAXDIM; i++) {
      arr[i][0] = 1;
      for (int j = 1; j <= i; j++) {
        arr[i][j] = arr[i-1][j-1] + arr[i-1][j];
      }
    }
  }
  return arr[n][r];
}

class Point {
public:
  struct badIndex { };
  struct dimMismatch { };

  Point(char const *commaSepCoords) {
    int *tmpcoord = new int[MAXDIM];
    nDim = 0;
    for (;;) {
      tmpcoord[nDim++] = atoi(commaSepCoords);
      if ((commaSepCoords = strchr(commaSepCoords, ',')) == 0) {
        break;
      }
      chkDim(nDim);
      commaSepCoords++;
    }
    coord = new int[nDim];
    memcpy(coord, tmpcoord, sizeof(int) * nDim);
  }

  Point(int dimArg) {
    nDim = dimArg;
    coord = new int[nDim];
  }

  Point(Point const& other, int index) {
    nDim = other.nDim;
    coord = new int[nDim];
    for (int i = nDim; --i >= 0;) {
      coord[i] = index % other[i];
      index /= other[i];
    }
  }

  Point(Point const& other) {
    nDim = other.nDim;
    coord = new int[nDim];
    for (int i = nDim; --i >= 0;) {
      coord[i] = other[i];
    }
  }

  void operator=(Point const& other) {
    if (this != &other) {
      delete coord;
      nDim = other.nDim;
      coord = new int[nDim];
      for (int i = nDim; --i >= 0;) {
	coord[i] = other[i];
      }
    }
  }

  int index(Point const& other) throw (dimMismatch) {
    if (other.nDim != nDim) {
      throw dimMismatch();
    }
    int retval = 0;
    for (int i = 0; i < nDim; i++) {
      retval = retval * other[i] + coord[i];
    }
    return retval;
  }

  ~Point() {
    delete coord;
  }
  int dim() const {
    return nDim;
  }
  int& operator[](int index) {
    if (index < 0 || index >= nDim) {
      throw badIndex();
    }
    return coord[index];
  }

  int const& operator[](int index) const {
    if (index < 0 || index >= nDim) {
      throw badIndex();
    }
    return coord[index];
  }

  void operator-=(Point const &other) throw (dimMismatch) {
    if (other.nDim != nDim) {
      throw dimMismatch();
    }
    for (int i = nDim; i-- > 0; ) {
      coord[i] -= other.coord[i];
    }
  }

  void operator+=(Point const &other) throw (dimMismatch) {
    if (other.nDim != nDim) {
      throw dimMismatch();
    }
    for (int i = nDim; i-- > 0; ) {
      coord[i] += other.coord[i];
    }
  }

  bool operator>(Point const &other) const throw (dimMismatch) {
    if (other.nDim != nDim) {
      throw dimMismatch();
    }
    for (int i = nDim; i-- > 0; ) {
      if (coord[i] <= other.coord[i]) {
        return false;
      }
    }
    return true;
  }
  bool operator>=(Point const &other) const throw (dimMismatch) {
    if (other.nDim != nDim) {
      throw dimMismatch();
    }
    for (int i = nDim; i-- > 0; ) {
      if (coord[i] < other.coord[i]) {
        return false;
      }
    }
    return true;
  }
  bool operator==(Point const &other) const throw (dimMismatch) {
    if (other.nDim != nDim) {
      throw dimMismatch();
    }
    for (int i = nDim; i-- > 0; ) {
      if (coord[i] != other.coord[i]) {
        return false;
      }
    }
    return true;
  }

  bool operator<(Point const &other) const throw (dimMismatch) {
    if (other.nDim != nDim) {
      throw dimMismatch();
    }
    for (int i = nDim; i-- > 0; ) {
      if (coord[i] >= other.coord[i]) {
        return false;
      }
    }
    return true;
  }

  bool operator<=(Point const &other) const throw (dimMismatch) {
    if (other.nDim != nDim) {
      throw dimMismatch();
    }
    for (int i = nDim; i-- > 0; ) {
      if (coord[i] > other.coord[i]) {
        return false;
      }
    }
    return true;
  }

  bool operator!=(Point const &other) const throw (dimMismatch) {
    return !(*this == other);
  }

  int taxi(Point const &other) const throw (dimMismatch) {
    if (other.nDim != nDim) {
      throw dimMismatch();
    }
    int retval = 0;
    for (int i = nDim; i-- > 0; ) {
      retval += abs(coord[i] - other.coord[i]);
    }
    return retval;
  }

  int taxi(Point const &maxCoords, int index) const throw (dimMismatch) {
    if (maxCoords.nDim != nDim) {
      throw dimMismatch();
    }
    int retval = 0;
    for (int i = nDim; --i >= 0;) {
      retval += abs(coord[i] - index % maxCoords[i]);
      index /= maxCoords[i];
    }
    return retval;
  }

private:
  int nDim;
  int *coord;
};

static ostream& operator<<(ostream& out, Point const &p) {
  for (int i = 0; i < p.dim(); i++) {
    if (i > 0) {
     out << ", ";
    }
    out << p[i];
  }
  return out;
}

//  This class solves a bunch of simultaneous equations by
//  Gauss-Jacobi elimination
//
class Gaussian {
public:
  Gaussian(int nvarArg) {
    nvar = nvarArg;
    int ncol = nvar + 1;
    storage = new double[nvar * ncol];
    rows = new double*[nvar];
    for (int i = 0; i < nvar; i++) {
      rows[i] = &storage[i * ncol];
    }
  }
  ~Gaussian() {
    delete storage;
    delete rows;
  }
  int size() const {
    return nvar;
  }
  double *operator[](int row) {
    return rows[row];
  }
  double const *operator[](int row) const {
    return rows[row];
  }
  void solve() {
    for (int row = 0; row < nvar; row++) {
      for (int col = nvar; col >= row; col--) {
        rows[row][col] /= rows[row][row];
      }
      for (int other = 0; other < nvar; other++) {
        if (other != row) {
	  double fac = rows[other][row];
	  for (int col = row; col <= nvar; col++) {
	    rows[other][col] -= rows[row][col] * fac;
	  }
	}
      }
      cerr << row << "...      \r";
    }
  }
private:
  int nvar;
  double **rows;
  double *storage;
};

static ostream& operator<<(ostream& out, Gaussian const &g) {
  for (int row = 0; row < g.size(); row++) {
    for (int col = 0; col <= g.size(); col++) {
      char buf[42];
      sprintf(buf, "%6.2f", g[row][col]);
      out << buf;
    }
    out << "\n";
  }
  return out;
}

//  This class represents a general resistance graph
//
class ResistGraph {
public:
  ResistGraph() {
    gaussian = 0;
  }
  virtual ~ResistGraph() {
    delete gaussian;
  }
  Gaussian *gaussian;
  virtual bool adjacent(int index1, int index2) const = 0;
  virtual bool init() {
    if (!internal_init()) {
      return false;
    }

    //  Create equation matrix
    gaussian = new Gaussian(nVert);
    Gaussian &g = *gaussian;
    for (int row = 0; row < nVert; row++) {
      if (row == test1Ind || row == test2Ind) {
	for (int col = 0; col < nVert; col++) {
	  g[row][col] = (row == col) ? 1 : 0;
	}
	g[row][nVert] = (row == test2Ind) ? 1 : 0;
      } else {
	g[row][row] = 0;
	for (int col = 0; col < nVert; col++) {
	  if (row != col) {
	    if (adjacent(row, col)) {
	      g[row][col] = -1;
	      g[row][row]++;
	    } else {
	      g[row][col] = 0;
	    }
	  }
	}
	g[row][nVert] = 0;
      }
    }

    return true;
  }
  int test1Ind;
  int test2Ind;
  int nVert;
  virtual bool theory(double &result) {
    return false;
  }
  double result() const {
    Gaussian &g = *gaussian;
    double res = 0;
    for (int col = 0; col < nVert; col++) {
      if (adjacent(test1Ind, col) == 1) {
	res += g[col][nVert];
      }
    }
    return 1/res;
  }
  virtual void graphDesc(ostream& out) = 0;
  virtual void printPoint(ostream& out, int index) = 0;
  virtual void testDesc(ostream& out) { }
  virtual void desc(ostream& out) {
    graphDesc(out);
    out << "  Test points at ";
    printPoint(out, test1Ind);
    out << " and ";
    printPoint(out, test2Ind);
    testDesc(out);
    out << "\n";
  }
protected:
  virtual bool internal_init() = 0;
};

class Lattice : public ResistGraph {
public:
  Lattice() :
    minP("0"), maxP("0"), test1P("0"), test2P("0"), dims("0")
  { }
  void setBounds (Point const &minPArg, Point const& maxPArg)
  {
    minP = minPArg;
    maxP = maxPArg;
    dims = maxP;
    dims -= minP;
    nVert = 1;
    for (int i = 0; i < dims.dim(); i++) {
      dims[i]++;
      nVert *= dims[i];
    }
  }
  void setTests(Point const &test1PArg, Point const &test2PArg) {
    test1P = test1PArg;
    test2P = test2PArg;
  }
  virtual void graphDesc(ostream& out) {
    out << "Lattice ";
    for (int i = 0; i < dims.dim(); i++) {
      if (i > 0) {
        out << "x";
      }
      out << dims[i];
    }
    out << " from " << minP << " to " << maxP << "\n";
  }
  virtual void printPoint(ostream& out, int index) {
    Point p(dims, index);
    p += minP;
    out << "(" << p << ")";
  }
  virtual bool adjacent(int index1, int index2) const {
    return Point(dims, index1).taxi(dims, index2) == 1;
  }
  virtual bool internal_init() {
    try {
      if (!(maxP >= minP)) {
	fmtError("Expected max-point coords to be >= min-point coords");
	return false;
      }
      if (!(test1P >= minP)) {
	fmtError("Expected test-point-1 coords to be >= min-point coords");
	return false;
      }
      if (!(test2P >= minP)) {
	fmtError("Expected test-point-2 coords to be >= min-point coords");
	return false;
      }
      if (!(test1P <= maxP)) {
	fmtError("Expected test-point-1 coords to be <= max-point coords");
	return false;
      }
      if (!(test2P <= maxP)) {
	fmtError("Expected test-point-2 coords to be <= max-point coords");
	return false;
      }
      if (test2P == test1P) {
	fmtError("Expected distinct test-points");
	return false;
      }
    } catch (Point::dimMismatch) {
      fmtError("Expected points to have same dimensions");
      return false;
    }

    test1P -= minP;
    test2P -= minP;

    test1Ind = test1P.index(dims);
    test2Ind = test2P.index(dims);

    test1P += minP;
    test2P += minP;

    return true;
  }
  Point minP;
  Point maxP;
  Point test1P;
  Point test2P;
  Point dims;
};

class Hypercube : public Lattice {
public:
  int dim;
  double theory_result;
  Hypercube(int dimArg) : dim(dimArg) {
    Point minPArg(dim);
    Point maxPArg(dim);
    for (int i = 0; i < dim; i++) {
      minPArg[i] = 0;
      maxPArg[i] = 1;
    }
    setBounds(minPArg, maxPArg);
    theory_result = 0;
  }
  void testOpp() {
    setTests(minP, maxP);
    theory_result = 0;
    for (int i = 0; i < dim; i++) {
      theory_result += 1.0/choose(dim - 1, i);
    }
    theory_result /= dim;
  }
  void testAdj() {
    Point adj(minP);
    adj[0] = 1;
    setTests(minP, adj);
    theory_result = ((double)((1 << dim) - 1))/(dim * (1 << (dim-1)));
  }
  virtual void graphDesc(ostream& out) {
    out << "Hypercube - " << dim << " dimensions\n";
  }
  virtual bool theory(double &result) {
    result = theory_result;
    return true;
  }
  virtual void testDesc(ostream& out) {
    out << " [";
    out << (adjacent(test1Ind, test2Ind) ? "adjacent" : "opposing");
    out << " vertices]";
  }
};

class Trigrid : public Lattice {
public:
  virtual void graphDesc(ostream& out) {
    out << "Triangular Grid ";
    for (int i = 0; i < dims.dim(); i++) {
      if (i > 0) {
        out << "x";
      }
      out << dims[i];
    }
    out << " from " << minP << " to " << maxP << "\n";
  }
  virtual bool adjacent(int index1, int index2) const {
    Point p1(dims, index1);
    Point p2(dims, index2);
    return p1.taxi(p2) == 1 ||
           (p1[0] == p2[0] - 1 && p1[1] == p2[1] + 1) ||
           (p1[0] == p2[0] + 1 && p1[1] == p2[1] - 1);
  }
  virtual bool internal_init() {
    if (!Lattice::internal_init()) {
      return false;
    }
    if (dims.dim() != 2) {
      fmtError("Expected 2D coords for TRIGRID");
      return false;
    }
    return true;
  }
};

class Hexgrid : public Lattice {
public:
  virtual void graphDesc(ostream& out) {
    out << "Hexagonal Grid ";
    for (int i = 0; i < dims.dim(); i++) {
      if (i > 0) {
        out << "x";
      }
      out << dims[i];
    }
    out << " from " << minP << " to " << maxP << "\n";
  }
  virtual bool adjacent(int index1, int index2) const {
    Point p1(dims, index1);
    Point p2(dims, index2);
    if (p1.taxi(p2) != 1) {
      return false;
    }
    if (p1[1] != p2[1]) {
      //  One unit apart, y coord different
      return true;
    }
    return (p1[0] < p2[0])
      ? ((p2[0]+p2[1]-minP[0]-minP[1]) & 1) != 0
      : ((p1[0]+p1[1]-minP[0]-minP[1]) & 1) != 0;
  }
  virtual bool internal_init() {
    if (!Lattice::internal_init()) {
      return false;
    }
    if (dims.dim() != 2) {
      fmtError("Expected 2D coords for HEXGRID");
      return false;
    }
    return true;
  }
};

class Simplex : public ResistGraph {
public:
  Simplex(int dimArg) : dim(dimArg) {
    test1Ind = 0;
    test2Ind = 1;
    nVert = dim + 1;
  }
  virtual bool adjacent(int index1, int index2) const {
    return (index1 == index2) ? 0 : 1;
  }
  virtual void graphDesc(ostream& out) {
    out << "Simplex - " << dim << " dimensions\n";
  }
  virtual void printPoint(ostream& out, int index) {
    out << "(";
    for (int i = 0; i < dim; i++) {
      int coord = (index == 0) ? 0 : (i == index - 1) ? 1 : 0;
      if (i > 0) {
	out << ", ";
      }
      out << coord;
    }
    out << ")";
  }
  int dim;
  virtual bool internal_init() {
    return true;
  }
  virtual bool theory(double &result) {
    result = 2.0 / (dim + 1);
    return true;
  }
};

class Crosspoly : public ResistGraph {
public:
  Crosspoly(int dimArg) : dim(dimArg) {
    nVert = 2 * dim;
  }
  virtual bool internal_init() {
    return true;
  }
  virtual bool adjacent(int index1, int index2) const {
    return (index1/2 == index2/2) ? 0 : 1;
  }
  virtual void graphDesc(ostream& out) {
    out << "Cross Polytope - " << dim << " dimensions\n";
  }
  virtual void printPoint(ostream& out, int index) {
    int axis = index / 2;
    bool neg = (index % 2 == 0);
    out << "(";
    for (int i = 0; i < dim; i++) {
      int coord = (i == axis) ? (neg ? -1 : 1) : 0;
      if (i > 0) {
	out << ", ";
      }
      out << coord;
    }
    out << ")";
  }
  void testOpp() {
    test1Ind = 0;
    test2Ind = 1;
  }
  void testAdj() {
    test1Ind = 0;
    test2Ind = 2;
  }
  virtual void testDesc(ostream& out) {
    out << " [";
    out << (adjacent(test1Ind, test2Ind) ? "adjacent" : "opposing");
    out << " vertices]";
  }
  virtual bool theory(double &result) {
    if (adjacent(test1Ind, test2Ind)) {
      result = (dim - 0.5) / (dim * (dim - 1));
    } else {
      result = 1.0/(dim - 1);
    }
    return true;
  }
  int dim;
};

int main(int argc, char *argv[]) {

  //  Set up program name
  //
  if ((progName = strchr(argv[0], '/')) == 0 &&
      (progName = strchr(argv[0], '\\')) == 0) {
    progName = argv[0];
  } else {
    progName++;
  }

  bool verbose = false;

  //  Handle command-line options
  //
  int currArg = 1;
  while (currArg < argc && argv[currArg][0] == '-') {
    switch (argv[currArg++][1]) {
    case 'v':
      verbose = true;
      break;
    default:
      usageError();
    }
  }

  //  Open input file
  FILE *fpIn;

  if (currArg == argc) {
    inFile = "<stdin>";
    fpIn = stdin;
  } else if (currArg + 1 == argc) {
    inFile = argv[currArg];
    if ((fpIn = fopen(inFile, "r")) == 0) {
      cerr << progName << ": Could not open file '"
           << inFile << "' for reading\n";
      exit(1);
    }
  } else {
    usageError();
  }

  char *whiteSpace = " \t\r\n";

  //  Process command lines
  //
  for (; fgets(lineBuf, sizeof(lineBuf), fpIn) != 0; lineNo++) {
    char *s = lineBuf;

    //  Skip initial blanks on the line
    //
    while (isspace(*s)) {
      s++;
    }

    //  Skip empty / comment lines
    //
    if (*s == '\0' || *s == '#') {
      continue;
    }

    //  Read keyword at beginning of line
    //
    char tmp[sizeof(lineBuf)];
    strcpy(tmp, s);
    s = tmp;
    if (!isalpha(*s)) {
      fmtError("Expected an alphabetic keyword");
      continue;
    }
    while (isalpha(*s)) {
      s++;
    }
    if (!isspace(*s)) {
      fmtError("Expected an alphabetic keyword followed by a space");
      continue;
    }
    *s++ = '\0';

    ResistGraph *graph = 0;

    if (strcasecmp(tmp, "lattice") == 0 ||
        strcasecmp(tmp, "trigrid") == 0 ||
        strcasecmp(tmp, "hexgrid") == 0) {
      char *rangelo = strtok(s, ";");
      char *rangehi = strtok(0, ";");
      char *test1 = strtok(0, ";");
      char *test2 = strtok(0, ";");
      char *leftover = strtok(0, ";");

      if (test2 == 0 || leftover != 0) {
	fmtError("Expected 3 semicolons on a data line");
	continue;
      }

      Lattice *lat;
      switch (tolower(tmp[0])) {
      case 't': lat = new Trigrid(); break;
      case 'h': lat = new Hexgrid(); break;
      default:  lat = new Lattice(); break;
      }
      lat->setBounds(rangelo, rangehi);
      lat->setTests(test1, test2);
      graph = lat;
    } else if (strcasecmp(tmp, "simplex") == 0) {
      char *dims = strtok(s, whiteSpace);
      char *leftover = strtok(0, whiteSpace);
      int ndim;
      if (dims == 0 || leftover != 0) {
	fmtError("Expected  dimensions  after keyword SIMPLEX");
	continue;
      } else if ((ndim = atoi(dims)) < 1) {
	fmtError("Expected dimension >=1 for simplex");
	continue;
      }
      chkDim(ndim);
      graph = new Simplex(ndim);
    } else if (strcasecmp(tmp, "hypercube") == 0) {
      char *dims = strtok(s, whiteSpace);
      char *oppadj = strtok(0, whiteSpace);
      char *leftover = strtok(0, whiteSpace);
      int ndim;
      if (oppadj == 0 || leftover != 0) {
	fmtError("Expected  dimensions  and  OPP|ADJ  after keyword HYPERCUBE");
	continue;
      } else if ((ndim = atoi(dims)) < 1) {
	fmtError("Expected dimension >=1 for hypercube");
	continue;
      }
      bool opp = false;
      if (strncasecmp(oppadj, "OPP", 3) == 0) {
        opp = true;
      } else if (strncasecmp(oppadj, "ADJ", 3) != 0) {
        fmtError("Expected trailing keyword 'OPP' or 'ADJ'");
	continue;
      }
      chkDim(ndim);
      Hypercube *hc = new Hypercube(ndim);
      if (opp) {
        hc->testOpp();
      } else {
        hc->testAdj();
      }
      graph = hc;
    } else if (strcasecmp(tmp, "crosspoly") == 0) {
      char *dims = strtok(s, whiteSpace);
      char *oppadj = strtok(0, whiteSpace);
      char *leftover = strtok(0, whiteSpace);
      int ndim;
      if (oppadj == 0 || leftover != 0) {
	fmtError("Expected  dimensions  and  OPP|ADJ  after keyword CROSSPOLY");
	continue;
      } else if ((ndim = atoi(dims)) < 2) {
	fmtError("Expected dimension >=2 for hypercube");
	continue;
      }
      bool opp = false;
      if (strncasecmp(oppadj, "OPP", 3) == 0) {
        opp = true;
      } else if (strncasecmp(oppadj, "ADJ", 3) != 0) {
        fmtError("Expected trailing keyword 'OPP' or 'ADJ'");
	continue;
      }
      chkDim(ndim);
      Crosspoly *cp = new Crosspoly(ndim);
      if (opp) {
        cp->testOpp();
      } else {
        cp->testAdj();
      }
      graph = cp;
    } else {
	fmtError("Unrecognised keyword");
	continue;
    }

    if (graph != 0) {
      if (graph->init()) {

	cerr << "solving (" << graph->nVert << " points) ...\n";
	graph->gaussian->solve();

	char result[20];
	sprintf(result, "%.10f", graph->result());
	graph->desc(cout);
	cout << "    Result by solving network = " << result << "\n";
	double theory;
	if (graph->theory(theory)) {
	  sprintf(result, "%.10f", theory);
	  cout << "    Result by theory          = " << result << "\n";
	}
	if (verbose) {
	  Gaussian &g = *graph->gaussian;
	  for (int i = 0; i < g.size(); i++) {
	    cerr << "    V";
	    graph->printPoint(cerr, i);
	    cerr << " = ";
	    sprintf(result, "%.10f", g[i][g.size()]);
	    cerr << result << "\n";
	  }
	}
      }
      delete graph;
    }
  }
  return 0;
}